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Abstract 

We analyze a large system of globally coupled phase oscillators whose natural frequencies are 
bimodally distributed. The dynamics of this system has been the subject of longstanding inter- 
est. In 1984 Kuramoto proposed several conjectures about its behavior; ten years later, Crawford 
obtained the first analytical results by means of a local center manifold calculation. Nevertheless, 
many questions have remained open, especially about the possibility of global bifurcations. Here 
we derive the system's stability diagram for the special case where the bimodal distribution consists 
of two equally weighted Lorentzians. Using an ansatz recently discovered by Ott and Antonsen, 
we show that in this case the infinite-dimensional problem reduces exactly to a flow in four dimen- 
sions. Depending on the parameters and initial conditions, the long-term dynamics evolves to one 
of three states: incoherence, where all the oscillators are desynchronized; partial synchrony, where 
a macroscopic group of phase-locked oscillators coexists with a sea of desynchronized ones; and a 
standing wave state, where two counter-rotating groups of phase-locked oscillators emerge. Ana- 
lytical results are presented for the bifurcation boundaries between these states. Similar results are 
also obtained for the case in which the bimodal distribution is given by the sum of two Gaussians. 

PACS numbers: 05.45.Xt 
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I. INTRODUCTION 



A. Background 

Large systems consisting of many coupled oscillatory units occur in a wide variety of 
situations Thus the study of the behaviors that such systems exhibit has been an 
active and continuing area of research. An important early contribution in this field was the 
introduction in 1975 by Kuramoto 3] of a simple model which illustrates striking features 
of such systems. Kuramoto employed two key simplifications in arriving at his model: (i) 
the coupling between units was chosen to be homogeneous and all-to-all (i.e., 'global'), so 
that each oscillator would have an equal effect on all other oscillators; and (ii) the oscillator 
states were solely described by a phase angle 6(t), so that their uncoupled dynamics obeyed 
the simple equation dOi/dt = uJi, where Ui is the intrinsic natural frequency of oscillator i, 

^ 1 is the number of oscillators, and i = 1,2, . . . , N. The natural frequencies Ui are, 
in general, different for each oscillator and are assumed to be drawn from some prescribed 
distribution function g^u). 

Much of the research on the Kuramoto mode 
unimodal (for reviews of this literature, see 
to be symmetric about a maximum at frequency u = ujq and to decrease monotonically 
and continuously to zero as \u! — u!q\ increases. In that case, it was found that as the 
coupling strength K between the oscillators increases from zero in the large- limit, there is 
a continuous transition at a critical coupling strength Kc = 2/{TTg{ujQ)). For K below K^, the 
average macroscopic, time-asymptotic behavior of the system is such that the oscillators in 
the system behave incoherently with respect to each other, and an order parameter (defined 
in Sec. [Tll) is correspondingly zero. As K increases past Kc, the oscillators begin to influence 
each other in such a way that there is collective global organization in the phases of the 
oscillators, and the time-asymptotic orderparameter assumes a non-zero constant value 
that increases continuously for K > Kc 

It is natural to ask how these results change if other forms of g{uj) are considered. In this 
paper we will address this question for what is perhaps the simplest choice of a non-unimodal 
frequency distribution: we consider a distribution g{uj) that has two peaks j^, 9| and is the 
sum of two identical unimodal distributions g, such that g{uj) = ^[g{uj — uj()) +g{Cd + LOQ)]. We 



las focused on the case where q(u) is 
1 

6j). Specifically, g is usually assumed 
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find that this modification to the original problem introduces qualitatively new behaviors. 
As might be expected, this problem has been previously addressed jsi 1^. However, due to 
its difficulty, the problem was not fully solved, and, as we shall show, notable features of the 
behavior were missed. 



B. Reduction method 



The development that makes our analysis possible is the recent paper of Ott and Anton- 



sen 



Using the method proposed in Ref. 



we reduce the original problem formulation 
from an integro-partial-differential equation [4, |5|, Q] for the oscillator distribution function (a 
function of cu, 9 and t) to a system of just a few ordinary differential equations (ODEs). Fur- 
thermore, we analyze the reduced ODE system to obtain its attractors and the bifurcations 
they experience with variation of system parameters. 

The reduced ODE system, however, represents a special restricted class of all the possible 
solutions of the original full system [Hj]. Thus a concern is that the reduced system might 
miss some of the actual system behavior. In order to check this, we have done numerical 
solutions of the full system. The result is that, in all cases tested, the time-asymptotic 
attracting behavior of the full system and the observed attractor bifurcations are all con- 
tained in, and are quantitatively described by, our ODE formulation. Indeed a similar result 
applies for the application of the method of Ref. to the original Kuramoto model with 
unimodally distributed frequencies 2], 3] and to the problem of the forced Kuramoto model 
with periodic drive 
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proven to capture 



12l |. Throughout 



On the other hand, the reduction method has not been mathematica 
all the attractors, for any of the systems to which it has been applied 
this paper we operate under the assumption (based on our numerical evidence) that the 
reduction method is reliable for the bimodal Kuramoto model. But we caution the reader 
that in general the situation is likely to be subtle and system-dependent; see Sec. IVID II for 
further discussion of the scope and limits of the reduction method. 
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C. Outline of the paper 



The organization of this paper is as follows. In Sec. |TT1 we formulate the problem and 
reduce it to the above-mentioned ODE description for the case where g{uj) is a sum of 
Cauchy-Lorentz distributions. 

Sec. mil provides an analysis of the ODE system. The main results of Sec. IIIII are a 
delineation of the different types of attractors that can exist, the regions of parameter space 
that they occupy (including the possibility of bistability and hysteresis), and the types of 
bifurcations that the attractors undergo. 

In Sec. \IV\ we establish that the attractors of the ODEs obtained in Section IIIII under 
certain symmetry assumptions are attractors of the full ODE system. In Section |Vl we 
confirm that these attractors and bifurcations are also present in the original system. In 
addition, we investigate the case where g{uj) is a sum of Gaussians, rather than Cauchy- 
Lorentz distributions. We find that the attractors and bifurcations in the Lorentzian case 
and in the Gaussian case are of the same types and that parameter space maps of the 
different behaviors are qualitatively similar for the two distributions. 

Finally, in Sec. |Vl] we compare our results to the earlier work of Kuramoto and 
Crawford [10|]. Then we discuss the scope and limits of the reduction method used here, and 



offer suggestions for future research. 

II. GOVERNING EQUATIONS 
A. Problem definition 

We study the Kuramoto problem of oscillators with natural frequencies Ui, 

^ =^. + f Esin(^,(t)-0„(t)), (1) 

where 6i are the phases of each individual oscillator and K is the coupling strength. We 
study this system in the limit N ^ oo for the case in which the distribution of natural 
frequencies is given by the sum of two Lorentzian distributions: 

A / 1 1 \ 

^^"^^ ~ 2^ V - ^o)' + A2 + (a; + u;o)2 + A2 J ■ 
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Here A is the width parameter (half-width at half-maximum) of each Lorentzian and icjo 
are their center frequencies, as displayed in Fig. [H A more physically relevant interpretation 
of ujq is as the detuning in the system (proportional to the separation between the two center 
frequencies). Note that we have written the distribution g{uj) so that it is symmetric about 



g(^) 



/' A 

/' A 

/' A 

/' A 

/' A 
/' A 
/' A 
/' * \ 

/' * V 


1 / \\ 

/' A 
/' A 
/' A 
/ ' A 
/' A 
/ ' A 
/ ' A 


/'2 A' \ 

// \ 


/ 2 a' \ 


/' \ 
/' ^ 


/ \ 

' \ 
\ \X 



FIG. 1: A bimodal distribution of natural frequencies, 5(0;), consisting of the sum of two 
Lorentzians. 

zero; this can be achieved without loss of generality by going into a suitable rotating frame. 

Another point to observe is that g{uj) is bimodal if and only if the peaks are sufficiently 
far apart compared to their widths. Specifically, one needs ujq > A/y/S. Otherwise the 
distribution is unimodal and the classical results of 



5[ would still apply. 



B. Derivation 

In the limit where N —* 00, Eq. ([1]) can be written in a continuous formulation m 
terms of a probability density f{6,uj,t). Here / is defined such that at time t, the fraction of 
oscillators with phases between 6 and 9 + dO and natural frequencies between uo and uj + duj 
is given by f{6, uj, t) dO duj. Thus 

'■2-K 

f{e,uo,t)deduo = 1 (3) 



00 Jo 



and 

r.27r 



f{e,u,t)de = g{u:), (4) 
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by definition of g{uj). 

The evolution of / is given by the continuity equation describing tlie conservation of 
oscillators: 



0, (5) 
where v{9,uj,t) is the angular velocity of the oscillators. From Eq. ([1]), we have 

(6) 



v{e,uj,t) = uj + K i f{e',uj,t)sm{e' ~ 

Jo 



Following Kuramoto, we define a complex order parameter 



Z{t) 



-oo JO 



e''f{e,Lu,t) dedtu 



(7) 



whose magnitude \z{t)\ < 1 characterizes the degree to which the oscillators are bunched 
in phase, and arg (z) describes the average phase angle of the oscillators. Expressing the 
velocity ([6]) in terms of z we obtain 



v{9,LJ,t) = u + Klm[ze 



(8) 
(9) 



where the * denotes complex co njug ate. 

Following Ott and Antonsen HI], we now restrict attention to a special class of density 
functions. By substituting a Fourier series of the form 
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l + 5^(/n(^,t)e^"' + c.c.) 



n=l 



(10) 



where 'c.c' stands for the complex conjugate of the preceeding term, and imposing the ansatz 
that 

Ucu,t)=a{u,tr, (11) 

we obtain 

da K 



where 



ot 2 



a{t, uj)g{uj)duj. 



(12) 



(13) 



We now consider solutions of (fT2l) and (fT3!l for initial conditions a{uj, 0) that satisfy the 
following additional conditions: (i) 10(^7,^)1 < 1; (ii) a{uj,0) is analytically continuable into 
the lower half plane lm{uj) < 0; and (iii) ^ as lm{uj) — oo. If these conditions 

are satisfied for a{uj,0), then, as shown in [11], they continue to be satisfied by a{u!,t) as it 
evolves under Eqs. (1121) and (IT^ . Expanding g{uj) in partial fractions as 

1 1 1 1 " 

{u — ujq) — iA {u — ujq) + iA {u + uq) — iA {u + ujq) + iAj' 

we find it has four simple poles at = icuo ± iA. Evaluating ( fT3l) by deforming the 
integration path from the real cu-axis to lm{uj) — oo, the order parameter becomes 

z{t) = liz,it) + Z2{t)), (14) 

where 

2i,2(t) = a*{±uJo-iA,t). (15) 

Substitution of this expression into (fT2l) yields two coupled complex ODEs, describing 
the evolution of two 'sub'-order parameters, 

zi = -(A + iuJo)zi 

+^[zi + Z2-{zl + z*)zl] (16) 



Z2 = -(A - iuJo)z2 

+ J [Zl + Z2-{ZI+Z*)zl] 



(17) 



where we use dots to represent the time derivative from now on. (This system agrees with the 
results of HI] for the case of two equal groups of oscillators with uniform coupling strength 
and average frequencies Uq and —ujq.) 



C. Reductions of the system 

The system derived so far is four-dimensional. If we introduce polar coordinates Zj = 
PjC^'^^ and define the phase difference ip = (p2 — (pi, the dimensionality can be reduced to 
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three: 



pi = -Api + ^(l-p?)(pi+P2Cos^) (18) 
K 

P2 = -Ap2 + — (1 - P2)(PlC0S^ + P2) (19) 

• K pI + pI + 2pjpl . 

= zcjo sm?/;. (20) 

4 P1P2 



To facilitate our analysis we now look for solutions of Eqs. fll8tl20p that satisfy the sym- 
metry condition 

pi(t) =P2(t) = Pit). (21) 

In Sec. HV] we will verify that these symmetric solutions are stable to perturbations away 
from the symmetry manifold and that the attractors of Eqs. (fT6| [T7|) lie within this manifold. 
Our analysis of the problem thus reduces to a study in the phase plane: 

P = (^1 - ^ - + (1 - P') cos^j (22) 
K 

^ = 2cJo - Y (1 + p^)sinV^. (23) 

III. BIFURCATION ANALYSIS 

Figure [2] summarizes the results of our analysis of Eqs. fl22l [23]) . We find that three types 



tnree types 



of attractors occur: the well-known incoherent and partially synchronized states 
corresponding to fixed points of (l22 l[23l) . as well as a standing wave state jlo| corresponding 
to limit-cycle solutions. In addition, we will show that the transitions between these states 
are mediated by transcritical, saddle-node, Hopf, and homoclinic bifurcations, as well as by 
three points of higher codimension. 



A. Scaling 

To ease the notation we begin by scaling Eqs. (122| [23ll . If we define q = p^ and non- 

dimensionalize the parameters and time such that 

K 

A ^ ^ (24, 

4u;o 

= ^ 
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FIG. 2: The bifurcation diagram for the Kuramoto system with a bimodal frequency distri- 
bution consisting of two equally weighted Lorentzians. The various bifurcation curves are de- 
noted as follows: TC=transcritical, SN=saddle-node, HB= (degenerate) Hopf, HC=homoclinic, 
and SNIPER=Saddle-node-infinite-period. The insets, labeled (a)-(g), show {q,ip) phase portraits 
in polar coordinates corresponding to the regions where the insets are located (see arrows for the 
boxed insets). Solid red dots and loops denote stable fixed point and limit cycles, respectively; open 
green dots are saddle fixed points, and open gray circles are repelling fixed points. All parameters 
refer to their original (unsealed) versions. 
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we obtain the dimensionless system 



ip = ljq ~ {1 + q) sinip. 



q = g(l — A — g + — g) cos ip) 



(25) 
(26) 



Here the overdot now means differentiation with respect to dimensionless time, and we 
have dropped all the tildes for convenience. For the rest of this section, all parameters 
will be assumed to be dimensionless (so there are implicitly tildes over them) unless stated 
otherwise. 

B. Bifurcations of the incoherent state 

The incoherent state is defined by pi = p2 = 0, or by g = in the phase plane formulation. 
The linearization of the incoherent state, however, is most easily performed in Cartesian 
coordinates using the formulation in Eqs. ( |T6l) and ( |T71) . We find the degenerate eigenvalues 



This degeneracy is expected because the origin is always a fixed point and because of the 
rotational invariance of that state. It follows that the incoherent state is stable if and only 
if the real parts of the eigenvalues are less than or equal to zero. 

The boundary of stable incoherence therefore occurs when the following conditions are 
met: 



These equations define the semicircle and the half-line shown in Fig. [2l labeled TC (for 
transcritical) and HB (for Hopf bifurcation), respectively. (Independent confirmation of 
these results can be obtained from the continuous formulation of Eq. ([1]) directly, as shown 
in the Appendix.) More precisely, we find that crossing the semicircle corresponds to a 
degenerate transcritical bifurcation, while crossing the half-line corresponds to a degenerate 
supercritical Hopf bifurcation. 

In the latter case, the associated limit-cycle oscillation indicates that the angle increases 
without bound; this reflects an increasing difference between the phases of the two 'sub'- 
order parameters of Eqs. (fT6l [T7|) . In terms of the original model, this means that the 




(27) 



(28) 




Q for uq < I 
for Uq > 1. 
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oscillator population splits into two counter-rotating groups, each consisting of a macroscopic 
number of oscillators with natural frequencies close to one of the two peaks of g{uj). Within 
each group the oscillators are frequency-locked. Outside the groups the oscillators remain 
desynchronized, drifting relative to one another and to the locked groups. This is the state 
™Q called a. wave. Intuitively speaking, it occurs when the two humps in 
the frequency distribution are sufficiently far apart relative to their widths. In Kuramoto's 
vivid terminology [sl, the population has spontaneously condensed into "a coupled pair of 
giant oscillators." 



C. Fixed point solutions and saddle- node bifurcations 

Along with the trivial incoherent state q = 0, the other fixed points of Eqs. fl25| [26]) 
satisfy 1 — A — g=(g — 1) cosip, and uq = {q + 1) sinip. Using trigonometric identities, we 
obtain 



ujq \^ f 1 — A — q^ 



q-1 



(29) 



or equivalently. 



1+q 



ooo = ±^3^VA(2-2g- A). (30) 



Thus, the fixed point surface q = q{ujQ, A) is defined implicitly. It can be single- or double- 
valued as a function of uiq for fixed A. To see this, consider how uq behaves as g ^ 0+. We 
find that 



VA(2 - A) 



(31) 



from which we observe that the behavior changes qualitatively at A = 3/2, as shown in 
Fig. [31 

The surface defined by p = p{ujq, A) can be plotted parametrically using p and A, as is 
seen in Fig. SJ The fold in the surface corresponds to a saddle-node bifurcation. Plots of the 
phase portrait of (g, tp) reveal that the upper branch of the double- valued surface in Fig. [3] 
corresponds to sinks, and the lower branch to saddle points; see Fig. [2] (c), (d), and (g). 

In physical terms, the sink represents a stable partially synchronized state, which is famil- 
iar from the classic Kuramoto model with a unimodal distribution 

mm 

. The oscillators 
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FIG. 3: Saddle-node bifurcation: at A = 3/2, q becomes double- valued. 

whose natural frequencies are closest to the center of the frequency distribution g{uj) become 
rigidly locked, and maintain constant phase relationships among themselves — in this sense, 
they act collectively like a "single giant oscillator," as Kuramoto [3;] put it. Meanwhile the 
oscillators in the tails of the distribution drift relative to the locked group, which is why one 
describes the synchronization as being only partial. 

The saddle points also represent partially synchronized states, though of course they 
are unstable. Nevertheless they play an important role in the dynamics because they can 
annihilate the stable partially synchronized states; this happens in a saddle-node bifurcation 
along the fold mentioned above. To calculate its location analytically, we use (l30l) and impose 
the condition for a turning point, duo/dq = 0, which yields 

g^-4g + 3-2A = 0. (32) 

Eliminating q from this equation using fl30l) . we obtain the equation for the saddle-node 
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FIG. 4: Fixed point surface. Bifurcation curves at the origin and the saddle-node curve are 
emphasized in black. 

bifurcation curve 

uo = ^J2 - lOA - A2 + 2(1 + 2A)3/2. (33) 

This curve is labeled SN in Fig.[2l Its intersection with the semicircle TC occurs at (cuq, A) = 
(^, |), and is labeled B in the figure. Note also that point C in the figure is not a Takens- 
Bagdanov point, as the saddle- node and Hopf bifurcations occur at different locations in the 
state space; see Figs. [2] (a) and (g). 

D. Bistability, homoclinic bifurcations, and SNIPER 

An examination of the dynamics corresponding to the approximately triangular param- 
eter space region ABC in Fig. [2] shows bistability. More specifically, we find that the stable 
incoherent fixed point coexists with the stable partially synchronized state produced by the 
saddle-node bifurcation described above, as shown in the state-space plot in Fig. ^c). 

Further study of these state-space plots led us to the homoclinic bifurcation curve marked 
HC, which was obtained numerically. The coexistence of states continues into region ACD, 
where we found that the stable partially synchronized state now coexists with the stable 
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limit cycle created at the Hopf curve. (See Fig. EJ^g).) This limit cycle is then destroyed by 
crossing the homoclinic curve, which is bounded by point A on one side and by point D on 
the other. 

At point D, the homoclinic curve merges with the saddle-node curve. This codimension- 
;wo bifurcation, occurring at approximately (1.3589, 0.7483), is known as a saddle- node- loop 
1^. Below D, however, the saddle-node curve exhibits an interesting feature: the saddle- 
node bifurcation occurs on an invariant closed curve. This bifurcation scenario is known as 
a saddle- node infinite-period bifurcation, or in short, SNIPER. If we traverse the SNIPER 
curve from left to right, the sink and saddle (the stable and unstable partially synchronized 
states) coalesce, creating a loop with infinite period. Beyond that, a stable limit cycle then 
appears — see Figs. [2] (d), (e), and (f). 

In conclusion, we have identified six distinct regions in parameter space and have identified 
the bifurcations that occur at the boundaries. 



IV. TRANSVERSE STABILITY 



Our analysis so far has been based on several simplifying assumptions. First, we restricted 
attention to a special family of oscillator distribution functions f{6,uj,t) and a bimodal 
Lorentzian form for g{uj), which enabled us to reduce the original infinite-dimensional system 
to a three-dimensional system of ODEs, Eqs. fll8ti20p . Second, we considered only symmetric 
solutions of these ODEs, by assuming pi = P2', this further decreased the dimensionality from 
three to two. 

The next two sections test the validity of these assumptions. We begin here by showing 
that the non-zero fixed point attractor (the stable partially synchronized state) and the limit 
cycle attractor (the standing wave state) for Eqns. f2E\} are transversely stable to small 
symmetry-breaking perturbations, i.e., perturbations off the invariant manifold defined by 
Pi = P2- This does not rule out the possible existence of attractors off this manifold, but it 
does mean that the attractors in the two-dimensional symmetric manifold are guaranteed 
to constitute attractors in the three-dimensional ODE system f fT8ll20l) . 

Let K = K/4 and consider the reduced governing equations ( fT8ll20l) without symmetry. 
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Introducing the longitudinal and transversal variables 

P± = ^(Pi-P2), (34) 

and substituting these into ( fT8ll20l) . we derive the equation for the transversal component 

P± = P± (k - A) - k{3pI + p]_) - Kcosip{l + pI - p]_] 

which describes the order parameter dynamics off the symmetric manifold. 

To simplify the notation, let q\\ = p| and q± = p]_ and scale the system using Eqs. ( l24ll . as 
before. Linearization and evaluation at the asymptotic solution denoted by (q'c^o); which 
may be either a fixed point or a limit cycle, yields the variational equation 

5q± = X±Sq± (35) 

where 

A_L = 1 - A - 3go - (1 + go) cosV'o- (36) 

Observe that 6q\i and Sip do not appear in linear order on the right hand side of fl35l) . This 
decoupling implies that A_l is the eigenvalue associated with the transverse perturbation 6q±, 
in the case where go is a fixed point. Similarly, if go is a limit cycle, the Floquet exponent 
associated with Sq± is simply (A_l), where the brackets denote a time average over one period. 
Hence the fixed point will be transversely stable if A_l < 0. The analogous condition for the 
limit cycle is (A_l) < 0. 

A. Fixed point stability 

To test the transverse stability of sinks for the two-dimensional flow, we solve Eq. (125|) 
for fixed points and obtain 

= 1 - A - go + (1 - go) cos^/^o- (37) 
Subtracting this from ( 136|) . we find 

A± = -2(go + cos?/^o)- (38) 
15 



Hence cosz/'o > is a sufficient condition for transverse stability. But at a non-trivial fixed 
point, 

cosV'o = ^ , (39) 

so the transverse stability condition is equivalent to go + ^ > 1- 

We claim that this inequality holds everywhere on the upper branch of the fixed point 
surface ( 130|) . Obviously the inequality is satisfied at all points where A > 1. For all other 



cases, consider the turning point from Fig. [3] defined by g^n = 2 ± vT+2A. Since the 
function of interest, Q(A) = qsn + A, has a global minimum with (5(0) = 1, and g^n is 
independent of ujq (at fixed A), it is a lower bound for all g(co'o) on the upper sheet of 
the fixed point surface, provided that g(co'o) is monotonically decreasing on the interval of 
[O,^^^^]. In fact, it is easier to establish that > duo/dq = A/D(g^ — 4g + 3 — 2A) with 
D = {q — 1)^a/2A — 2gA — A^; the latter expression is positive, and g^ — 4g + 3 — 2 A < 
whenever 1 > g > qsn- Thus transverse stability for the nodes on the fixed point surface 
follows. 

B. Limit cycle stability 

To examine the transverse linear stability of the limit cycle, we calculate the transverse 
Floquet exponent by averaging the eigenvalue over the period of one oscillation: 

(Ai) = 1 - A - 3 (go) - ((cosz^o) + (goCosV^o)) • (40) 

In order to render this expression definite, we rewrite Eq. fl25l) in terms of the limit cycle 
solution (go, ^ipo): 

^ (In go) = 1 - A - go + (1 - go) cos tpo. (41) 
Periodicity on the limit cycle guarantees (^Ingo) = 0, and so we have 

= 1 - A - (go) + ((1 - go) cos^^o), (42) 
which we subtract from the averaged eigenvalue to yield 

(Ax) = -2((go) + (cos^o)). (43) 
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Although we are not able to analytically demonstrate that (A_l) in (H3|l is negative, we have 
calculated (go) and (cos'?/'o) numerically for the limit cycle attractors of Eqs. (11811201) . This 
was done for 2500 parameter values corresponding to a grid in dimensionless parameter space, 
by sampling 50 evenly spaced values uj G [0.01,2.5] and A G [0.01,2.1]. The simulations 
were run with = 1024 oscillators. In all the cases that we tested, we found that (A_l) < 0. 

V. NUMERICAL EXPERIMENTS 

All of the results described above were obtained using the reduced ODE models derived 
in Sec. |TT] B and C, and are therefore subject to the restrictions described therein. It is 
therefore reasonable to ask if these results agree with the dynamics of the original system 
given in Eq. ([1]). To check this, a series of direct simulations of Eq. ([1]) using A^ = 10,000 
oscillators and fourth-order Runge-Kutta numerical integration were performed. 

First, we compared solutions of Eq. ([T]) with those of our reduced system Eqs. ( I22l [23|l 
in the region where we predicted the coexistence of attractors. For example, we show in 
Fig. [5] a bifurcation diagram computed along the line Auq/K = 1.092 that traverses the 
region ABCD in Fig. [2l (Note that here and for the rest of the paper, we revert to using 
the original, dimensional form of the variables.) The vertical lines in Fig. [5] indicate the 
locations of the bifurcations that were identified using the ODE models. For each point 
plotted, the simulation was run until the order parameter exhibited its time-asymptotic 
behavior; this was then averaged over the subsequent 5000 time steps. Error bars denote 
standard deviation. Note in particular the hysteresis, as well as the point with the large 
error bar, indicating the predicted limit cycle behavior. 

Next, we examined the behavior of Eq. ([1]) at 121 parameter values corresponding to an 
11 X 11 regular grid superimposed on Fig. [21 ranging from 0.1 to 2.1 at intervals of 0.2 on 
each axis. (In all cases, K was set to 1, and A and ujq were varied.) An additional series 
was run using a smaller grid (from 0.6 to 1.6 at intervals of 0.1 on each axis), to focus on 
the vicinity of region ABCD in Fig. [2l Initial conditions were chosen systematically in 13 
different ways, as follows: 

1. The oscillator phases were uniformly distributed around the circle, so that the overall 
order parameter had magnitude r = 0. 
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FIG. 5: Hysteresis loop as observed when traversing the bistable regions shown in Fig. [5] in the 
directions shown (arrows) along the line at Alvq/K = 1.092. The data were obtained from a 
simulation of Equation ([1]) with N = 10, 000 and K = 1. Vertical lines indicate where the reduced 
ODE models of Section [II] predict homoclinic (HC), degenerate Hopf (HB), and saddle-node (SN) 
bifurcations. Note that the point marked 'limit cycle' has a large error bar, reflecting the oscillations 
in the order parameter. 

2. The oscillators were all placed in phase at the same randomly chosen angle in [0, 27r], 
so that r = 1. 

3. The remaining 11 initial conditions were chosen by regarding the system as composed 
of two sub-populations, one for each Lorentzian in the bimodal distribution of frequen- 
cies, as in jsl. In one of the sub-populations, the initial phases of the oscillators were 
chosen to be randomly spaced within the angular sector [c + d,c — d]^ where c was 
chosen randomly in [0, 2tt] and d was chosen at random such that the sub-order pa- 
rameter magnitude ri = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, or 0.9 (all approximately). 
The result was that ri had one of these magnitudes and its phase was random in 
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[0, 27i]. The same procedure was followed for the other sub-population, subject to the 
constraint that ri 7^ r2. Our idea here was to deliberately break the symmetry of the 
system initially, to test whether it would be attracted back to the symmetric subspace 
defined by Eq. mh. 




FIG. 6: (a) The bifurcation diagram for the Kuramoto system with a bimodal frequency distribu- 
tion consisting of two equally weighted Gaussians. All the features in Fig. [2] are present, but are 
somewhat distorted. The transcritical (TC) and (degenerate) Hopf curves (HB) were obtained as 
described in the Appendix. The dotted lines represent conjectured saddle-node, homoclinic, and 
SNIPER curves. These are based on the numerically-observed bifurcations shown in (b), which is a 
magnification of the central region of (a). The symbols represent saddle-node (circles), homoclinic 
(triangles), and SNIPER (squares) bifurcations. 

In all the cases we examined, no discrepancies were found between the simulations and 
the predicted behavior. Although these tests were not exhaustive, and certainly do not 
constitute a mathematical proof, they are consistent with the conjecture that no additional 
attractors beyond those described in Section III exist. 

We then investigated the generality of our results by replacing the bimodal Lorentzian 
natural frequency distribution, Eq. ([2]), with the sum of two Gaussians: 

g{uj) = —= e + e (44) 
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and computing the corresponding bifurcation diagram analogous to Fig. [2l The results 
are shown in Fig. [61 The transcritical (TC) and degenerate Hopf bifurcation (HB) curves 
were obtained using the continuous formulation of Eq. ([T]); see the Appendix for details. 
In addition, saddle-node, homoclinic, and SNIPER bifurcations were numerically observed 
at several parameter values, and based on these data, we estimated the location of the 
corresponding curves (dashed lines). All the features of Fig. [2] are preserved, but the curves 
are somewhat distorted. 



VI. DISCUSSION 



We conclude by relating our work to three previous studies, and then offer suggestions 
for further research, both theoretical and experimental. 



A. Kuramoto's conjectures 

In his book on coupled oscillators, Kuramoto [s^ speculated about how the transition 
from incoherence to mutual synchronization might be modified if the oscillators' natural 
frequencies were bimodally distributed across the population. On pp. 75-76 of Ref. js], he 
wrote "So far, the nucleation has been supposed to be initiated at the center of symmetry 
of g. This does not seem to be true, however, when g is concave there." His reasoning was 
that for a bimodal system, synchrony would be more likely to start at the peaks of g. If 
that were true, it would mean that a system with two equal peaks would go directly from 
incoherence to having two synchronized clusters of oscillators, or what we have called the 
standing wave state, as the coupling K is increased. The critical coupling at which this 
transition would occur, he argued, should be Kc = 2/ {ng{LJms.x)), analogous to his earlier 
result for the unimodal case. According to this scenario, the synchronized clusters would 
be tiny at onset, comprised only of oscillators with natural frequencies near the peaks of 
g{uj). Because of their small size, Kuramoto claimed these clusters "will behave almost 
independently of each other." With further increases in K, however, the clusters "will come 
to behave like a coupled pair of giant oscillators, and for even stronger coupling they will 
eventually be entrained to each other to form a single giant oscillator." (This is what we 
have called the partially synchronized state.) 
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Let us now re-examine Kuramoto's conjectures in light of our analytical and numerical 
results, as summarized in Fig. [7|(a). For a fair comparison, we must assume that g is con- 
cave at its center frequency = 0; for the bimodal Lorentzian (Eq. ([2]), this is equivalent 
to ciJo/A > (Otherwise g is unimodal and incoherence bifurcates to partial synchro- 

nization as K is increased, consistent with Kuramoto's classic result as well as the lowest 
portion of Fig. [D^a).) 

So restricting attention from now on to the upper part of FigI7](a) where cuq/A > l/-\/3, 
what actually happens as K increases? Was Kuramoto right that the bifurcation sequence 
is always incoherence — >■ standing wave — partial synchronization? 

No. For tuo/A between 1/v^ and 1 (meaning the distribution is just barely bimodal), in- 
coherence bifurcates directly to partial synchronization — the "single giant oscillator" state — 
without ever passing through an intermediate standing wave state. In effect, the system still 
behaves as if it were unimodal. But there is one new wrinkle: we now see hysteresis in 
the transition between incoherence and partial synchronization, as reflected by the lower 
bistable region in Fig. [7](a). 

Is there any part of Fig. [Tl^a) where Kuramoto's scenario really does occur? Yes — but 
it requires that the peaks of g be sufficiently well separated. Specifically, suppose cJo/A > 
1.81 . . ., the value at the codimension-2 saddle-node-loop point where the homoclinic and 
SNIPER curves meet (i.e., point D in Fig. [2]). In this regime everything behaves as Kuramoto 
predicted. 

An additional subtlety occurs in the intermediate regime where the peaks of g are neither 
too far apart nor too close together. Suppose that 1 < c^o/A < 1.81 . . .. Here the system 
shows a different form of hysteresis. The bifurcations occur in the sequence that Kuramoto 
guessed as K increases, but not on the return path. Instead, the system skips the stand- 
ing wave state and dissolves directly from partial synchronization to incoherence as K is 
decreased. 

Finally we note that Kuramoto's conjectured formula Kc = 2/{'Kg{ujjn&x)) is incorrect, 
although it becomes asymptotically valid in the limit of widely separated peaks. Specifically, 
his prediction is equivalent to K^, = , ~ 4A(1 — jid^/ujof ), which approaches 

l+-\/l+(A/i^o)^ 

the correct result K^. = 4 A as ujq/A oo. 
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B. Crawford's center manifold analysis 



Crawford lOj obtained the first mathematical results for the system studied in this paper. 
Using center manifold theory, he calculated the weakly nonlinear behavior of the infinite- 
dimensional system in the neighborhood of the incoherent state. From this he derived the 
stability boundary of incoherence. His analysis also included the effects of white noise in 
the governing equations. 
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FIG. 7: Left: Results from our analysis, white: incoherence, dark gray: partial synchronization, 
light gray: standing wave (limit cycles), vertical lines: coexistence of incoherent and partially 
synchronized states, horizontal lines: coexistence of partial synchronization and standing waves. 
Right: Crawford's bifurcation diagram in [10(]. In our study there is no noise, and so the diffusion 
is = 0. Crawford's e corresponds to our A. /: Incoherent states, PS: partially synchronized, 
SW: standing wave, equivalent to what we describe as two counterrotating flocks of oscillators. 
(Permission to print by Springer Verlag.) 



Figure [7](b), reproduced from Fig. 4 in Ref. [lOj], summarizes Crawford's findings. Here 
D is the noise strength (note: our analysis is limited to D = 0), e is the width of the 
Lorentzians (equivalent to A in our notation), and icuo are the center frequencies of the 
Lorentzians (as here). The dashed line in Fig. [Tl^b) shows Crawford's schematic depiction of 
the unknown stability boundary between the standing waves and the partially synchronized 
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state. He suggested a strategy for calculating this boundary, and highlighted it as an open 
problem, writing in the figure caption, "...the precise nature and location of this boundary 
have not been determined." Our results, summarized in Figs. [2] and Fig. [7](a), now fill in the 
parts that were missing from Crawford's analysis. 



C. Stochastic model of Bonilla et al. 

In a series of papers (see jo] for a review), Bonilla and his colleagues have explored 
what happens if one replaces the Lorentzians in the frequency distribution with 5-functions, 
and adds white noise to the governing equations. The resulting system can be viewed as a 
stochastic counterpart of the model studied here; in effect, the noise blurs the (5-functions into 
bell-shaped distributions analogous to Lorentzians or Guassians. And indeed, the system 
shows much of the same phenomenology as seen here: incoherence, partially synchronized 
states, standing waves, and bistability [6|]. 

However, a complete bifurcation diagram analogous to Fig. [2] has not yet been worked out 
for this model. The difficulty is that no counterpart of the ansatz (ITT!) has been found; the 
stochastic problem is governed by a second-order Fokker-Planck equation, not a first-order 
continuity equation, and the Ott-Antonsen ansatz (fTTj) no longer works in this case. Perhaps 
there is some way to generalize the ansatz appropriately so as to reduce the stochastic model 
to a low-dimensional system, but for now this remains an open problem. 



D. Directions for future research 



There are several other questions suggested by the work described here. 



1. Validity of reduction method 

The most important open problem is to clarify the scope and limits of the Ott-Antonsen 
method used in Sec. Ill B[ Under what conditions is it valid to assume that the infinite- 
dimensional Kuramoto model can be replaced by the low-dimensional dynamical system 
implied by the Ott-Antonsen ansatz? Or to ask it another way, when do all the attractors of 
the infinite-dimensional system lie in the low-dimensional invariant manifold corresponding 
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to this ansatz? 

This question has now become particularly pressing, because two counterexamples have 
recently come to light in which the Ott-Antonsen method gives an incomplete account 
of the full system's dynamics. When the method was applied to the problem of chimera 
states for two interacting populations of identical phase oscillators, it predicted only sta- 
tionary and periodic chimeras [l^, whereas subsequent numerical experiments revealed that 



quasiperiodic chimeras can also exist and be stable 



15j. Likewise, chaotic states are known 



to emerge from a wide class of initial conditions for series arrays of identical overdamped 



Josephson junctions coupled through a resistive load [la, llTI]. Yet the Ott-Antonsen ansatz 



cannot account for these chaotic states, because the reduced ODE system turns out to be 



only two dimensional 
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19|. 



What makes this all the more puzzling is that the method works so well in other cases. 
It seems to give a full inventory of the attractors for the bimodal Kuramoto model studied 



here, as well as for the unimodal Kuramoto model in its original form 



external periodic forcing H 



ll| or with 



121, 



So we are left in the unsatisfying position of not knowing when the method works, or 
why. In some cases it (apparently) captures all the attractors, while in other cases it does 
not. How does one make sense of all this? 

A possible clue is that in all the cases where the method has so far been successful, the 
individual oscillators were chosen to have randomly distributed frequencies; whereas in the 
cases where it failed, the oscillators were identical. Perhaps the mixing induced by frequency 
dispersion is somehow relevant here? 

A resolution of these issues may come from a new analytical approach. Pikovsky and 



Rosenblum [15| and MiroUo, Marvel and Strogatz [18| have independently shown how to 
place the Ott-Antonsen ansatz [ll] in a more general mathematical framework by relating 
it to the group of Mobius transformations 18|, l20|] or, equivalently, to a trig onometric trans- 
formation ISj originally introduced in the study of Josephson arrays 17|. This approach 
includes the Ott-Antonsen ansatz as a special case, but is more powerful in the sense that 
it provably captures all the dynamics of the full system, and it works for any A^, not just 
in the infinite- limit. The drawback is that the analysis becomes more complicated. It 
remains to be seen what conclusions can be drawn — and, perhaps, what longstanding prob- 
lems can be solved — when this new approach is unleashed on the Kuramoto model and its 
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many relatives. 

Even in those instances where the Ott-Antonsen ansatz doesn't account for all the attrac- 
tors of the full system, it can still provide useful information, for instance by giving at least 
some of the attractors and by easing the calculation of them. Moreover, the transient evo- 
lution from initial conditions off the Ott-Antonsen invariant manifold can yield interesting 
phenomena not captured by the ansatz, as discussed in Appendix C of [21I. 



2. Asymmetric bimodal distributions 

Now returning to the specific problem of the bimodal Kuramoto model: What happens if 
the humps in the bimodal distribution have unequal weights? The analysis could proceed as 
in this paper, up to the point where we assumed symmetry between the two sub-populations. 
One would expect new phenomena such as traveling waves to arise because of the broken 
symmetry. 



3. Finite-size effects 



We have focused here exclusively on the infinite-A^ limit of the Kuramoto model. What 
happens when the number of oscillators is reduced? How do finite-size effects influence the 



bifurcation diagram? An analysis along the lines of 
these questions. 



23| could be fruitful for investigating 



4- Comparison with experiment 

Finally, it would be interesting to test some of these theoretical ideas in real systems. 
One promising candidate is the electrochemical oscillator system studied by Hudson and 



colleagues 2J], in which the frequency distribution can be bimodal or even multimodal 25|. 
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APPENDIX A: ALTERNATIVE CALCULATION OF THE BOUNDARY OF 
STABILITY FOR THE INCOHERENT STATE 



The system in Eq. ([T]), together with the bimodal natural frequency distribution given 
in Eq. ([2]), can be expressed using the formulation in 8[] as two interacting populations of 
oscillators. In this case, each population has a separate Lorenzian frequency distribution 
of width A and center frequency at Uo or —Uq, and the two- by-two matrix describing the 
relative coupling weights (i.e, Eq. (1) in j^) has 1/2 in each entry. By postulating that a 
small perturbation to the incoherent state grows exponentially as e'^*, and setting s = iv for 
the marginally stable state, Eq. (9) of Ref. [sf] gives the following expression for the critical 
coupling value K: 

^ ^ 2(A^-z.^+a;g)-Fz(4A^) 

The boundary of stability of the incoherent state is obtained by requiring that this expression 
be strictly real. One solution is obtained for z/ = 0, resulting in K = 2(A^ + ci;q)/A, which 
is equivalent to 

1 + -° =1. (A2) 



K J \K ^ 

This is the equation for the semicircle in Figure [21 corresponding to a transcritical bifurcation 
of the incoherent state. Another solution, obtained by assuming that z/ 7^ in Eq. flAip and 
requiring ujq> K = 4A. This is the equation for the half-line in Figure [2] corresponding 
to the degenerate Hopf bifurcation of the incoherent state. 

If the bimodal natural frequency distribution is given by a sum of Gaussians of standard 
deviation a and centers at icuo, then the two-population approach outlined above leads to 
the following equation: 

(A3) 



K = 

TT 




y2a 7 V 72' 



where 
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is known as the Faddeeva function and can be computed numerically 
requiring that K be real, two branches corresponding to v being equal and not equal to 
zero can be obtained. These are the boundaries of stability of the incoherent state shown in 
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